Load the necessary libraries
library(car) # for regression diagnostics
library(broom) # for tidy output
library(broom.mixed) ## for tidying mixed effects models
library(ggfortify) # for model diagnostics
library(sjPlot) # for outputs
library(knitr) # for kable
library(effects) # for partial effects plots
library(emmeans) # for estimating marginal means
library(ggeffects) # for partial effects plots
library(MASS) # for glm.nb
library(MuMIn) # for AICc
library(nlme)
library(lme4) # for lmer
library(lmerTest) # for satterthwaite p-values with lmer
library(performance) # for residuals diagnostics
library(see) # for plotting residuals
# library(pbkrtest) #for kenward-roger p-values with lmer
library(glmmTMB) # for glmmTMB
library(DHARMa) # for residuals and diagnostics
library(tidyverse) # for data wrangling
A plant pathologist wanted to examine the effects of two different strengths of tobacco virus on the number of lesions on tobacco leaves. She knew from pilot studies that leaves were inherently very variable in response to the virus. In an attempt to account for this leaf to leaf variability, both treatments were applied to each leaf. Eight individual leaves were divided in half, with half of each leaf inoculated with weak strength virus and the other half inoculated with strong virus. So the leaves were blocks and each treatment was represented once in each block. A completely randomised design would have had 16 leaves, with 8 whole leaves randomly allocated to each treatment.
Tobacco plant
Sampling design
Format of tobacco.csv data files
| LEAF | TREAT | NUMBER |
|---|---|---|
| 1 | Strong | 35.898 |
| 1 | Week | 25.02 |
| 2 | Strong | 34.118 |
| 2 | Week | 23.167 |
| 3 | Strong | 35.702 |
| 3 | Week | 24.122 |
| ... | ... | ... |
| LEAF | The blocking factor - Factor B |
| TREAT | Categorical representation of the strength of the tobacco virus - main factor of interest Factor A |
| NUMBER | Number of lesions on that part of the tobacco leaf - response variable |
tobacco <- read_csv("../public/data/tobacco.csv", trim_ws = TRUE)
glimpse(tobacco)
## Rows: 16
## Columns: 3
## $ LEAF <chr> "L1", "L1", "L2", "L2", "L3", "L3", "L4", "L4", "L5", "L5", …
## $ TREATMENT <chr> "Strong", "Weak", "Strong", "Weak", "Strong", "Weak", "Stron…
## $ NUMBER <dbl> 35.89776, 25.01984, 34.11786, 23.16740, 35.70215, 24.12191, …
The response here represents the average number of lesions across a number of microscope views within a section of leaf. As such, provided the averages are sufficiently far from 0, there is no a priori reason to suspect that the data would not follow a Gaussian distribution.
Model formula: \[ y_i \sim{} \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]
where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of the treatment on the number of lesions. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with leaves.
We will start by explicitly declaring the categorical variable (TREATMENT) as a factor. In addition, random effects (in this case LEAF) should also be declared as factors.
tobacco <- tobacco %>%
mutate(
LEAF = factor(LEAF),
TREATMENT = factor(TREATMENT)
)
To explore the assumptions of homogeneity of variance and normality, a boxplot of each Treatment level is appropriate.
ggplot(tobacco, aes(y = NUMBER, x = TREATMENT)) +
geom_boxplot()
Conclusions:
It can also be useful to get a sense of the consistency across blocks (LEAF). That is, do all Leaves have a similar baseline level of lesion susceptibility and do they respond similarly to the treatment.
ggplot(tobacco, aes(y = NUMBER, x = as.numeric(LEAF))) +
geom_line(aes(linetype = TREATMENT))
## If we want to retain the original LEAF labels
ggplot(tobacco, aes(y = NUMBER, x = as.numeric(LEAF))) +
geom_blank(aes(x = LEAF)) +
geom_line(aes(linetype = TREATMENT))
Conclusions:
Given that there are only two levels of Treatment (Strong and Weak), it might be easier to visualise the differences in baselines and effect consistency by plotting as:
ggplot(tobacco, aes(y = NUMBER, x = TREATMENT, group = LEAF)) +
geom_point() +
geom_line(aes(x = as.numeric(TREATMENT)))
Conclusions:
The above figure also serves as a good way to visualise certain aspects of mixed effects models. When we fit a mixed effects model that includes a random blocking effect (in this case LEAF), we are indicating that we are allowing there to be a different intercept for each block (LEAF). In the current case, the intercept will represent the first Treatment level (Strong). So the random effect is specifying that the intercept can vary from Leaf to Leaf.
We can think of the model as having two tiers (a hierarchy), where the tiers of the hierarchy represent progressively smaller (typically) spatial scales. In the current example, the largest spatial units are the leaves (blocking factor). Within the leaves, there are the two Treatments (Strong and Weak) and within the Treatments are the individual observations.
We tend to represent this hierarchy upside down in the model formula:
\[ y_i \sim{} \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i =\beta_0 + \boldsymbol{\beta} \bf{X_i}\\ \beta_0 = \boldsymbol{\gamma} \bf{Z_i} \]
In addition to allowing there to be a different intercept per leaf, we could allow there to be a different magnitude of effect (difference between Strong and Week Treatment) per leaf. That is considered a random slope. From the figure above, there is some evidence that the effects (slopes) may vary from Leaf to Leaf.
Incorporating a random slope (in addition to a random intercept), may reduce the amount of unexplained variance and thus improve the power of the main effect (Treatment effect).
As with Generalized Linear Models (GLM), Generalized Linear Mixed Models (GLMM) are fitted using maximum likelihood. However, for models containing random effects, maximum likelihood estimates are known to be biased. To correct this bias, rather than use maximum likelihood, we can use Residual Maximum Likelihood (REML). REML maximises the likelihood of the residuals instead of maximising the likelihood of the observed data.
## Fit the random intercepts model
tobacco.lme <- lme(NUMBER ~ TREATMENT, random = ~ 1 | LEAF, data = tobacco, method = "REML")
## Fit the random intercepts/slope model
tobacco.lme1 <- lme(NUMBER ~ TREATMENT, random = ~ TREATMENT | LEAF, data = tobacco, method = "REML")
Having now fit the two models (a random intercept model and a random intercept/slope model), we can use AIC (or AICc) to determine which model is ‘best’.
AICc(tobacco.lme, tobacco.lme1)
Conclusions:
Alternatively, we could use a likelihood ratio test (LRT).
anova(tobacco.lme, tobacco.lme1) %>% print()
## Model df AIC BIC logLik Test L.Ratio p-value
## tobacco.lme 1 4 96.70799 99.26422 -44.35400
## tobacco.lme1 2 6 100.67049 104.50484 -44.33525 1 vs 2 0.03749921 0.9814
Conclusions:
tobacco.lme - the random intercepts model).## Fit the random intercepts model
tobacco.lmer <- lmer(NUMBER ~ TREATMENT + (1 | LEAF), data = tobacco, REML = TRUE)
## Fit the random intercepts/slope model
## Note the following could not be run in lmer as there was not enough observations
## to estimate all of the effects and random effects.
## tobacco.lmer1 <- lmer(NUMBER ~ TREATMENT + (TREATMENT|LEAF), data=tobacco, REML=TRUE)
## tobacco.lmer1 <- lmer(NUMBER ~ TREATMENT + (TREATMENT|LEAF), data=tobacco, REML=TRUE,
## control = lmerControl(
## check.nobs.vs.nRE = "ignore",
## optimizer = 'optim',
## optCtrl = list(method = 'BFGS'))
## )
If we had been able to fit both models (the random intercept model and the random intercept/slope model), we could have used AIC (or AICc) to determine which model is ‘best’. As it is, we only have the random intercepts model.
# AICc(tobacco.lme, tobacco.lme1)
Similarly, we could have used a likelihood ratio test (LRT).
# anova(tobacco.lme, tobacco.lme1) %>% print
Optimizers include:
nlminb (non linear optimiser)Nelder-Mead - robust but relatively slowBFGS - quasi-Newton method (acronym of authors) - gradient based
L-BFGS-B - as above but with different constraintsSANN - simulated annealing - relatively slowglmmTMB uses the nlminb (nonlinear minimizer)
## Fit the random intercepts model
tobacco.glmmTMB <- glmmTMB(NUMBER ~ TREATMENT + (1 | LEAF),
data = tobacco, REML = TRUE
)
## Fit the random intercepts/slope model
## Note the following model did not converge - probably due to insufficient data.
tobacco.glmmTMB1 <- glmmTMB(NUMBER ~ TREATMENT + (TREATMENT | LEAF),
data = tobacco, REML = TRUE
)
## Try a different optimizer (BFGS)
tobacco.glmmTMB1 <- glmmTMB(NUMBER ~ TREATMENT + (TREATMENT | LEAF),
data = tobacco, REML = TRUE,
control = glmmTMBControl(
optimizer = "optim",
optArgs = "BFGS"
)
)
## optArgs = 'Nelder-Mead'))
If we had been able to fit both models (the random intercept model and the random intercept/slope model), we could have used AIC (or AICc) to determine which model is ‘best’. As it is, we only have the random intercepts model.
AICc(tobacco.glmmTMB, tobacco.glmmTMB1)
Similarly, we could have used a likelihood ratio test (LRT).
anova(tobacco.glmmTMB, tobacco.glmmTMB1) %>% print()
## Data: tobacco
## Models:
## tobacco.glmmTMB: NUMBER ~ TREATMENT + (1 | LEAF), zi=~0, disp=~1
## tobacco.glmmTMB1: NUMBER ~ TREATMENT + (TREATMENT | LEAF), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## tobacco.glmmTMB 4 96.708 99.798 -44.354 88.708
## tobacco.glmmTMB1 6 100.670 105.306 -44.335 88.670 0.0375 2 0.9814
# autoplot(tobacco.lme)
tobacco.lme %>%
plot_model(type = "diag") %>%
plot_grid()
Conclusions:
tobacco.lme %>% performance::check_model()
## Could not compute standard errors from random effects for diagnostic plot.
## Homogeneity of variance could not be computed. Cannot extract residual variance from objects of class 'lme'.
Conclusions:
## unfortunately, nlme is not supported...
# tobacco.resid = simulateResiduals(tobacco.lme, plot=TRUE)
plot_model() and plot_grid() are of limited value for mixed effects models since:
plot_grid()plot_model(tobacco.lmer, type = "diag")[-2] %>% plot_grid()
Conclusions:
tobacco.lmer %>% performance::check_model()
Conclusions:
tobacco.resid <- tobacco.lmer %>% simulateResiduals(plot = TRUE)
Conclusions:
plot_model() and plot_grid() are of limited value for mixed effects models since:
plot_grid()plot_model(tobacco.glmmTMB, type = "diag")[-2] %>% plot_grid()
Conclusions:
tobacco.glmmTMB %>% performance::check_model()
Conclusions:
tobacco.resid <- tobacco.glmmTMB %>% simulateResiduals(plot = TRUE)
Conclusions:
Unfortunately, plot_model does not support lme models.
# plot_model(tobacco.lme, type='eff')
tobacco.lme %>%
allEffects() %>%
plot()
There is a rather obtuse issue with ggpredict. Due to downstream processing (using predict()), models fit with data in the form of a tibble will not work. In order to use ggpredict on a lme model, it is necessary to ensure that the data are data.frame. We can simply update the model..
tobacco.lme1 <- update(tobacco.lme, data = as.data.frame(tobacco))
tobacco.lme1 %>%
ggpredict(terms = "TREATMENT") %>%
plot()
tobacco.lme1 <- update(tobacco.lme, data = as.data.frame(tobacco))
tobacco.lme1 %>%
ggemmeans(~TREATMENT) %>%
plot()
tobacco.lmer %>% plot_model(type = "eff")
## $TREATMENT
tobacco.lmer %>%
allEffects() %>%
plot()
tobacco.lmer %>%
ggpredict(terms = "TREATMENT") %>%
plot()
tobacco.lmer %>%
ggemmeans(~TREATMENT) %>%
plot()
tobacco.glmmTMB %>% plot_model(type = "eff")
## $TREATMENT
tobacco.glmmTMB %>%
allEffects() %>%
plot()
tobacco.glmmTMB %>%
ggpredict(terms = "TREATMENT") %>%
plot()
tobacco.glmmTMB %>%
ggemmeans(~TREATMENT) %>%
plot()
tobacco.lme %>% summary()
## Linear mixed-effects model fit by REML
## Data: tobacco
## AIC BIC logLik
## 96.70799 99.26422 -44.354
##
## Random effects:
## Formula: ~1 | LEAF
## (Intercept) Residual
## StdDev: 3.692213 3.802922
##
## Fixed effects: NUMBER ~ TREATMENT
## Value Std.Error DF t-value p-value
## (Intercept) 34.94013 1.873988 7 18.644797 0.0000
## TREATMENTWeak -7.87938 1.901461 7 -4.143856 0.0043
## Correlation:
## (Intr)
## TREATMENTWeak -0.507
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.61850365 -0.48452438 0.01133404 0.40899997 1.42778194
##
## Number of Observations: 16
## Number of Groups: 8
## to get confidence intervals
tobacco.lme %>% intervals()
## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## (Intercept) 30.50885 34.940130 39.37141
## TREATMENTWeak -12.37562 -7.879381 -3.38314
##
## Random Effects:
## Level: LEAF
## lower est. upper
## sd((Intercept)) 1.580532 3.692213 8.625218
##
## Within-group standard error:
## lower est. upper
## 2.252287 3.802922 6.421125
Conclusions:
tobacco.lme %>% tidy(effects = "fixed", conf.int = TRUE)
## including the random effects
tobacco.lme %>% tidy(conf.int = TRUE)
Conclusions:
# warning this is only appropriate for html output
tobacco.lme %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)
| Â | NUMBER | |||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| (Intercept) | 34.94 | 1.87 | 30.51 – 39.37 | <0.001 |
| TREATMENT [Weak] | -7.88 | 1.90 | -12.38 – -3.38 | 0.004 |
| Random Effects | ||||
| σ2 | 14.46 | |||
| τ00 LEAF | 13.63 | |||
| N LEAF | 8 | |||
| Observations | 16 | |||
| Marginal R2 / Conditional R2 | 0.534 / NA | |||
| AIC | 96.708 | |||
tobacco.lme %>% r.squaredGLMM()
## R2m R2c
## [1,] 0.3707883 0.6761022
## Nakagawa's R2
tobacco.lme %>% performance::r2_nakagawa()
## Random effect variances not available. Returned R2 does not account for random effects.
## # R2 for Mixed Models
##
## Conditional R2: NA
## Marginal R2: 0.534
Conclusions:
tobacco.lmer %>% summary()
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: NUMBER ~ TREATMENT + (1 | LEAF)
## Data: tobacco
##
## REML criterion at convergence: 88.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.61850 -0.48453 0.01133 0.40900 1.42778
##
## Random effects:
## Groups Name Variance Std.Dev.
## LEAF (Intercept) 13.63 3.692
## Residual 14.46 3.803
## Number of obs: 16, groups: LEAF, 8
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 34.940 1.874 11.332 18.645 7.38e-10 ***
## TREATMENTWeak -7.879 1.901 7.000 -4.144 0.00433 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TREATMENTWk -0.507
## to get confidence intervals
tobacco.lmer %>% confint()
## 2.5 % 97.5 %
## .sig01 0.000000 7.258832
## .sigma 2.333872 6.361451
## (Intercept) 31.216405 38.663855
## TREATMENTWeak -11.829135 -3.929629
Conclusions:
tobacco.lmer %>% tidy(effects = "fixed", conf.int = TRUE)
## including the random effects
tobacco.lmer %>% tidy(conf.int = TRUE)
Conclusions:
# warning this is only appropriate for html output
tobacco.lmer %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)
| Â | NUMBER | |||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| (Intercept) | 34.94 | 1.87 | 30.86 – 39.02 | <0.001 |
| TREATMENT [Weak] | -7.88 | 1.90 | -12.02 – -3.74 | 0.001 |
| Random Effects | ||||
| σ2 | 14.46 | |||
| τ00 LEAF | 13.63 | |||
| ICC | 0.49 | |||
| N LEAF | 8 | |||
| Observations | 16 | |||
| Marginal R2 / Conditional R2 | 0.371 / 0.676 | |||
| AIC | 96.708 | |||
tobacco.lmer %>% r.squaredGLMM()
## R2m R2c
## [1,] 0.3707882 0.6761026
## Nakagawa's R2
tobacco.glmmTMB %>% performance::r2_nakagawa()
## # R2 for Mixed Models
##
## Conditional R2: 0.676
## Marginal R2: 0.371
Conclusions:
tobacco.glmmTMB %>% summary()
## Family: gaussian ( identity )
## Formula: NUMBER ~ TREATMENT + (1 | LEAF)
## Data: tobacco
##
## AIC BIC logLik deviance df.resid
## 96.7 99.8 -44.4 88.7 14
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## LEAF (Intercept) 13.63 3.692
## Residual 14.46 3.803
## Number of obs: 16, groups: LEAF, 8
##
## Dispersion estimate for gaussian family (sigma^2): 14.5
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 34.940 1.874 18.645 < 2e-16 ***
## TREATMENTWeak -7.879 1.901 -4.144 3.42e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tobacco.glmmTMB %>% vcov()
## Conditional model:
## (Intercept) TREATMENTWeak
## (Intercept) 3.511833 -1.807776
## TREATMENTWeak -1.807776 3.615552
## the following is not compatible with piping
cov2cor(vcov(tobacco.glmmTMB)$cond)
## (Intercept) TREATMENTWeak
## (Intercept) 1.0000000 -0.5073298
## TREATMENTWeak -0.5073298 1.0000000
## to get confidence intervals
tobacco.glmmTMB %>% confint()
## 2.5 % 97.5 % Estimate
## cond.(Intercept) 31.267180 38.613080 34.940130
## cond.TREATMENTWeak -11.606175 -4.152587 -7.879381
## LEAF.cond.Std.Dev.(Intercept) 1.580539 8.625194 3.692216
tobacco.glmmTMB %>% ranef()
## $LEAF
## (Intercept)
## L1 -0.3539126
## L2 -1.5406161
## L3 -0.7111775
## L4 -2.5604485
## L5 -1.2399299
## L6 3.5023104
## L7 5.6216063
## L8 -2.7178321
Conclusions:
tobacco.glmmTMB %>% tidy(effects = "fixed", conf.int = TRUE)
## including the random effects
tobacco.glmmTMB %>% tidy(conf.int = TRUE)
Conclusions:
# warning this is only appropriate for html output
tobacco.glmmTMB %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)
| Â | NUMBER | |||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| (Intercept) | 34.94 | 1.87 | 31.27 – 38.61 | <0.001 |
| TREATMENT [Weak] | -7.88 | 1.90 | -11.61 – -4.15 | <0.001 |
| Random Effects | ||||
| σ2 | 14.46 | |||
| τ00 LEAF | 13.63 | |||
| ICC | 0.49 | |||
| N LEAF | 8 | |||
| Observations | 16 | |||
| Marginal R2 / Conditional R2 | 0.371 / 0.676 | |||
| AIC | 96.708 | |||
tobacco.glmmTMB %>% r.squaredGLMM()
## R2m R2c
## [1,] 0.3707882 0.6761025
## Nakagawa's R2
tobacco.glmmTMB %>% performance::r2_nakagawa()
## # R2 for Mixed Models
##
## Conditional R2: 0.676
## Marginal R2: 0.371
Conclusions:
Given that there was only two levels of a single predictor, the only follow up analysis might be to explore hypotheses around effects of specific magnitudes. For example, we could investigate whether the strong inoculation treatment increases the number of lesions by over a specific amount. For this example, lets explore the likelihood that the Strong treatment could result in an increase in lesions of over 2 (on average) compared to the Weak treatment.
tobacco.lme %>%
multcomp::glht(linfct = c("TREATMENTWeak=-2")) %>%
summary()
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lme.formula(fixed = NUMBER ~ TREATMENT, data = tobacco, random = ~1 |
## LEAF, method = "REML")
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## TREATMENTWeak == -2 -7.879 1.901 -3.092 0.00199 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
Conclusions:
## lme not supported
# brms::hypothesis(tobacco.lme, 'TREATMENTWeak < -2')
tobacco.lmer %>%
multcomp::glht(linfct = c("TREATMENTWeak=-2")) %>%
summary()
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lmer(formula = NUMBER ~ TREATMENT + (1 | LEAF), data = tobacco,
## REML = TRUE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## TREATMENTWeak == -2 -7.879 1.901 -3.092 0.00199 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
Conclusions:
tobacco.lmer %>% brms::hypothesis("TREATMENTWeak < -2")
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (TREATMENTWeak)-(-2) < 0 -5.91 1.87 -8.9 -2.87 499
## Post.Prob Star
## 1 1 *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
Conclusions:
tobacco.glmmTMB %>%
multcomp::glht(linfct = c("TREATMENTWeak=-2")) %>%
summary()
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glmmTMB(formula = NUMBER ~ TREATMENT + (1 | LEAF), data = tobacco,
## REML = TRUE, ziformula = ~0, dispformula = ~1)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## TREATMENTWeak == -2 -7.879 1.901 -3.092 0.00199 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
Conclusions:
## Does not appear to work for glmmTMB?
# brms::hypothesis(tobacco.glmmTMB, 'TREATMENTWeak < -2')
tobacco.lme %>%
emmeans(~TREATMENT) %>%
as.data.frame() %>%
ggplot() +
geom_pointrange(aes(y = emmean, x = TREATMENT, ymin = lower.CL, ymax = upper.CL))
tobacco.lmer %>%
emmeans(~TREATMENT) %>%
as.data.frame() %>%
ggplot() +
geom_pointrange(aes(y = emmean, x = TREATMENT, ymin = lower.CL, ymax = upper.CL))
tobacco.glmmTMB %>%
emmeans(~TREATMENT) %>%
as.data.frame() %>%
ggplot() +
geom_pointrange(aes(y = emmean, x = TREATMENT, ymin = lower.CL, ymax = upper.CL))